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We propose a new method to investigate the thermal properties of QCD with a small quark 
chemical potential ^. Derivatives of quark and gluonic observables with respect to /i are computed 
at = for 2 flavors of p-4 improved staggered fermions with ma — 0.1,0.2 on a 16^ x 4 lattice, 
and used to calculate the leading order Taylor expansion in /i of the location of the pseudocritical 
04 ■ point about = 0. This expansion should be well-behaved for the small values of /iq/^c ~ 0.1 

' relevant for RHIC phenomenology, and predicts a critical curve T'c(m) in reasonable agreement with 

, estimates obtained using exact reweighting. In addition, we contrast the case of isoscalar and 

isovector chemical potentials, quantify the effect of 7^ on the equation of state, and comment on 
^ 1' the complex phase of the fermion determinant in QCD with /i 7^ 0. 
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I. INTRODUCTION 



The study of the phase structure of QCD at non-zero temperature and baryon density is one of the most interesting 
topics in contemporary physics. Heavy- ion collision experiments are running at BNL and CERN with the goal of the 
experimental production of a new state of matter, the quark-gluon plasma On the theoretical side, novel color 
, superconducting and superfluid phases have been conjectured at high baryon densities For these reasons the need 
for numerical studies of the QCD phase transition using lattice gauge theory simulations, currently the most powerful 
^ quantitative approach to QCD, with both temperature T 7^ and quark chemical potential /Zq 7^ 0, is more urgent 
' than ever. Precise theoretical inputs from simulations in the vicinity of the QCD phase transition are indispensable 
I to the understanding of the heavy-ion collision experiments. 

Over the last several years, the numerical study of lattice QCD has been successful at zero chemical potential and 
high temperature j^. In contrast, because the quark determinant is complex at 7^ and Monte-Carlo simulation 
not directly applicable, studies at non-zero /i are still largely exploratory. Recent developments with /i 7^ can be 
classified in two categories [Q. At the low temperatures and high densities where the new phases are expected, studies 
of model field theories such as two-color (SU(2)) QCD and the NJL model have been made. The simulation is possible 
; ^ ■ because in both cases the quark determinant is positive definite so that conventional Monte Carlo methods can be 
\ used. The other case is high temperature and low density, which is phenomenologically more important for RHIC, 
since the QCD phase transition both in the early universe and in the interesting regime for heavy-ion collisions is 
expected at rather low density, e.g., /j,q ~ 15MeV (/iq/Tc ~ 0.1) for RHIC In this region the reweighting method, 
in which observables at /i 7^ are computed by performing simulations at Re(/i) = 0, is applicable Using this 
method, the first results on the phase structure in the (/i,T) plane were recently obtained by Fodor and Katz 0. 
Unfortunately, although in principle with infinite statistics this method is exact, rather general arguments suggest that 
in practise the region of applicability of the reweighting method becomes narrower as the lattice volume is increased. 
Another efficient method at low density is via a Taylor expansion obtained by computing the derivatives of physical 
quantities with respect to /i at = 0. This approach is not restricted to small lattices, because it requires only the 
expectation values of local fermion bilinears; these are measured effectively on large systems using stochastic methods, 
and might even be expected to self-average as the volume increases. Since analyticity is required, however, the values 
of /i which can be reached must be bounded by, e.g. the critical point expected in the (/i, T) plane for QCD with 2 
light flavors. Pioneering work in such a framework have been made by developing expansions for free energy, yielding 
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quark number susceptibility , for hadronic screening masses and in the context of the 3-dimensional effective 
theory 0. ^ ^ 

In this study, we investigate the transition temperature Tc as a function of ^ ^ 0. In Sec. |^ we propose a new 
method to compute derivatives of physical quantities with respect to /i. Details of our simulations performed on a 
16'^ X 4 lattice with quark masses m — 0.1,0.2 are presented in Sec. |l^. In Sec. ^ we check the feasibility of the 
method by calculating the derivative of the transition point with respect to m. Our main result, the calculation of the 
second derivative of /3c with respect to ^ for 2 flavor QCD, is given in Sec. 0. Using data on the lattice beta-function, 
we are then able to translate this result into physical units, yielding an estimate for the phase transition line Tc(/i). 
We also discuss the response of the pressure p{T) and energy density e{T) to non-zero /i, and estimate their variation 
along the critical line. Finally in this section we discuss the problem of the complex phase of the quark determinant. 



and show that the sign problem is mild in the region of the phase diagram relevant for RHIC physics. Section VI 
presents our conclusions. 



II. REWEIGHTING METHOD FOR THE ^-DIRECTION 

Ferrenberg and Swendsen's reweighting method is a very useful technique to investigate critical phenomena fl^ . In 
QCD the expectation value of an observable 0(/3, m,/i) can in principle be computed by simulation at (/3o,too,^o) 
using the following identity: 

= Z(/3,m,^) / 2'^C)(det Af(™,M))"^'e-^«(« (1) 

/Q^aNiiln dctM(m,p)-ln dot M(mo ,Aio)) g-Sg (/3)+Sg (/3o) \ 



^gaATfClndot M(m,/j)-lndct M (ma ,fia)) q-S g{l3)+Sg{f3o)^ 



(/3o,T?Xo,/Jo) 



Here M is the quark matrix, Sg the gauge action, iVf the number of flavors, and a = 1 or 1/4 for Wilson or staggered 
lattice fermions respectively. The chemical potential parameter /x — /iqO, where a is the lattice spacing. Because 
det M(/i) is complex for Re(/x) ^ 0, the expectation values in (Q) can only be estimated by conventional Monte Carlo 
importance sampling if the simulation is performed for /io zero or pure imaginary. Most of the attempts to calculate 
at /i 7^ have used variants of this method The reweighting factor for the gauge part is easy to compute by 
measuring the plaquette P^^, since 

-Sg{(3) + SgiP,) = (fi- f3o) P^-^^)^ (3) 

for the standard Wilson action, and extensions for improved actions are easy to derive. However, to compute the 
fermion part, the calculation of the fermion determinant is required for each point (m, /i) we want to study. Such a 
calculation is quite expensive and difficult to perform in practice. Fodor and Katz have performed such calculations, 
and by reweighting in both /i and (3 have succeeded in tracing out the critical line PdfJ-) and locating the critical 
endpoint on small lattices 0. Their method exploits the fact that the overlap between ensembles at different points 
along the coexistence line separating hadronic and quark-gluon plasma phases remains reasonably large on finite 
systems. 

Another problem of the reweighting method is the sign problem, which will be discussed in detail later. As /i increases 
from zero, the calculation of eqn.(^ becomes more difficult due to fluctuations in the phase of the denominator. To 
avoid these problems, we restrict ourselves to calculating derivatives of physical quantities with respect to /i, which 
can be done at fi — 0. This yields estimates of the physical quantity as a continuous function of fi in a narrow range 
of /i, but the region of applicability is not restricted to the immediate neighbourhood of the phase transition. This 
permits the development of a Taylor expansion of observables in powers of /i = fJ.qa; strictly speaking, in fact, the 
physically relevant expansion parameter which ultimately must govern convergence is the fugacity fiq/T = Ntfi. The 
Taylor expansion for the fermionic part of the reweighting factor around /i = is 

, /detM(fi)\ ^ /z" 9"lndetAf(0) „ 
aNf In ^aNty ^ ^ = > 7^„u". (4) 

^ ^ ' ^ n— 1 ^ n— 1 

We similarly expand fermionic observables such as the chiral condensate, 

(^V) = {^l X Nt)-^aN,{iTM~^), (5) 
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where the lattice size is x Nt, once again obtaining a continuous function for small /i. Using the formula 

expressions for 9"(lndet Af)/9/i" and 9"(trA/^^)/9/i" in terms of traces over products of local operators and inverse 
matrices can be developed: 

aindetM _ / _i9A/\ 
a^lndctM / ^ i^^^/fx / aM,^ 



93lndetM , a^MX ^ i9M_ , ^^^/X „ / ^ i , ^ i , ^ , aA/\ 



-tr M-^-—M-^ + 2tr M-^-—M~^-—M^^ 



+3tr Af"^— — Af-i— ^A//"i - 6tr Ar^ —-M-^ —-Nr^ —-M~^ . (8) 



- 6tr M-^ ^—M-^ ^—M-^ ^—M~^ 

\ Ofl Ofl J 



We apply the random noise method to calculate the derivatives of In det M and tr A^ ^ , which enables us to compute 
on rather large volumes in comparison with usual studies of QCD with /i ^ 0. Using sets of random noise vectors 
rjai which satisfy the condition: limjv„^oo(l/-^n) vliVaj = Sij, we rewrite the trace of products of dM/dfi and 

Ai'-i as 

tr Af-i- Af-M = lim —Yvl^ Af-^- M'^ria. 9 

M^^rja = X and M^^{dM/d^) ■ ■ - rja = x are obtained by solving AIx = rja or AIx — (dM/dfi) ■ ■ - rja, and we compute 
the RHS of eqn.(^ with finite TVn. The error for estimates of physical observables made from iVconf configurations 
is expected to decrease as (iVniVconf)"^^^- Further notes on the application of the noise method are given in the 
Appendix. 

By using the derivatives of both the reweighting factor and fermionic observable up to n-th order in /x, we can 
obtain the correct answer for the expectation value up to n-th order, which can be easily checked by performing a 
Taylor expansion of the expectation value, eqn.(|l|), directly for each physical observable. Of course, for a pure gluonic 
observable such as the Polyakov loop L only the expansion of Indet Af is needed. Furthermore, we should note that 
a,t n — the odd order derivatives of both In det A/ and trM~^ are pure imaginary and the even order derivatives are 
real. This property is proved using the identities for the fermion matrix: 

d" /Vf 1 d" A/f 

AftO,)=r5Af(-M)r5, and ^-_(^) ^ (-I)"r5^(-M)r5, (10) 

where Fs is 75 for Wilson fermions and (_i)a:i+2;2+2;3+2;4 Jq^. staggered. Then, at /i = 

tr ( M~'%-^M-^%-^M~^ ■■■] = (-l)"i+"^+-tr ( M'^^-^M'^^-^M'^ ■ 
Because the terms in the n-th derivative satisfy ni + n2 + • • • = n, we obtain 



(11) 



a"lndetA/\* _ „a"lndetAf 

di? J ' dj? ' ^ ^ 



dfi^ J ' ' dn" 
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Using this property and the fact that Z is a real function of /3, m and /x, we can exphcitly confirm that, if the operator 
has the property such that even order derivatives are real and odd order derivatives are pure imaginary at ^ = 0, 
e.g., iiJip) or its susceptibility, then all odd order derivatives of the expectation value of a physical quantity are zero 
at /.t = 0, as we expect from the symmetry under changing /i to — /i. The derivative of the expectation value can 
be written as a sum of products of expectation values composed of the operator, the reweighting factor and their 
derivatives, and the total number of differentiations in each term has to be odd for an odd order derivative. Hence 
all terms for odd derivatives contain at least one expectation value of a pure imaginary operator and hence vanish, 
since the expectation value of a pure imaginary operator is zero. Therefore the first non-trivial order of corrections 
to e.g. {ipil>) or its susceptbility, which we compute in this study is 0{^^); the truncation errors, so far unquantified, 
are 0{p^). 



Lilt: xcx_yiui t;Apa.iioiuii ui ciii upcicLtui uy ^y^— q ^nA^ • tu w^/-'"^ 

expression (0) for {0)(_p^i^i) can be rewritten 



In order to be more specific, let us define the Taylor expansion of an operator by J2^=o ^"M"- Then to 0{fi^) the 



^ ((00 + Oifi + exp(7^l/i + 7^2A^^ - ASg)) 

^ '^^^^^^ (exp(7^l^i + 7^2M2-A5,)) ' 

where expectation values on the RHS are measured with respect to an ensemble generated at (/3o , 0) . Extension of 
this formula to combining data from several ensembles using multi- histogramming is straightforward [l^ . Further 
details on the evaluation of (^4|) using the noise method for fermionic operators are given in the Appendix. 
In order to determine the pseudocritical point, we calculate the Polyakov loop susceptibility, 

XL=N^s{{L^)-{L?), (15) 
where the Polyakov loop L = (N^)^^ J2x -^c^^'^Ut Ui{x,t), and the susceptibihty of the chiral condensate!, 

X^^ = (iV3 X N,)-\aN,f (((trM-i)2) - {irM-'f) . (16) 
We define the transition point (3c{ii) by the peak position of these susceptibilities for each fi; 



d(3 



= 0. (17) 



If we compute dx/df) correctly up to n-th order in /j,, we can determine the n-th derivative of (3c with respect to 
/I. For example if we determine (3c{fJ.) using an operator such as (iptp)-, which is real and whose first derivative at 
mu=0 is pure imaginary, then the first derivative Pc{fJ.) vanishes because as argued above the first derivative of the 
susceptibility is zero in this case. 

Finally, note we can also estimate the magnitude of fluctuations of the phase of det M, because on each configuration 
this phase can be expressed in terms of the odd terms of the Taylor expansion of In det M; this will be discussed in 



more detail in section V C 



III. SIMULATIONS FOR iVp = 2 IMPROVED STAGGERED FERMIONS 

We e mploy a combination of the Symanzik improved gauge and 2 flavors of the p4-improved staggered fermion 



actions [b_3 14 1. The partition function is deflned by 



Z(/3,m,/z)= I VU{AetM)^'/^e-^\ (18) 
M. " 



y X, fl>L^ X, fl,!^ ) 



'^Note that we only calculate the disconnected part of the complete chiral susceptibility. 
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i+2j,y ^ij * '^j)^x-i-2lv 



+i-2j,y ^ij 



+2j,y 



+1+24, a 



-2prr(l,2) 



i~2i,y 



i+24,j/ 



E '5.+4+2^,, - e-^f/i;,^^) V - 4 - 2i)5^_i_^^y 

i 



(20) 



where W^^^(x) and W^^'^{x) are 1x1 and 1x2 Wilson loops, 77^(2;) = (— 1)^^^ jg the KS phase and 

f/^iV^a^) = ^ [U^.{x)Ul{x + fL-i')Ul[x + fi- 2v) + C/t(x - z>)C/t(x - 2z>)C/^(x - 2z))] , (21) 



1 



l + 6w 



C/^(a;) + w E [t^^(2;)C^A.(a; + i>)?7j;(a; + A) + Vl{x - v)U ^,{x - />)[/^(a; + /*-!>)] 



The coefScients are /? = 6/17^, ci = —1/12, cq = 1 — 8ci, cf = 3/8, cf = 1/96, and lo = 0.2. The action is derived 
such that rotational invariance of the free fermion propagator is restored up to 0(ji^\ It is known that this action 
makes the discretization error of the equation of state pressure p(T) small as T — ^ cx3, and obtained by this action 
is consistent with that obtained using improved Wilson fermions ||l^,^. To incorporate the chemical potential, we 
generalise the standard prescription of treating [i as an imaginary gauge potential Aq p^ ] by multipUjing the terms 
generating n-step hops in the positive and negative temporal directions by e"'^ and e""'' respectively □. 

We investigated the transition points for quark masses m = 0.1 and 0.2. The corresponding pseudoscalar and vector 
meson mass ratios are mps/mv ~ 0.70 and 0.85 [|4j. We compute the Polyakov loop, chiral condensate, and their 
susceptibilities. The simulations were performed on a 16'^ x 4 lattice for 7 values of /3 G [3.64,3.67] for m = 0.1 and 
6 values of /3 G [3.74, 3.80] for m = 0.2, using the Hybrid R algorithm. We adopted a step size At = 0.25 x m and a 
molecular dynamics trajectory length t = 0.5. For each trajectory 10 sets of Z2 noise vectors were used to calculate 
the reweighting factor and the derivatives of ipif) up to second order in /i. 

For the calculation of mass reweighting surveyed in Sec. we took a total of 220600 trajectories at m = 0.1 and 
155000 trajectories at m = 0.2. For the study with /i 7^ described in Sec. we used 128000 trajectories at to = 0.1 
and 86000 trajectories at m = 0.2. The details are summarized in Table ||. The multi-histogram method of [|l^ was 
used to reweight in the /3-direction using data from several values of (3. Errors were estimated using the jack-knife 
method with bin size 100 trajectories. 

IV. REWEIGHTING FOR QUARK MASS 

Before calculating derivatives with respect to /i, it is worthwhile to calculate the derivatives with respect to quark 
mass TO, which is not only potentially important for the chiral extrapolation, but also a good demonstration of the 
reweighting technique for a parameter appearing in the fermion action. Because we cannot compare the result obtained 
by reweighting in the n direction with the result of an actual simulation at /i 7^ 0, this test is a necessary check of 
the reliability of our method. The reweighting formula for quark mass is easily obtained from eqn. (^) and eqns. (|^) 



■^Note that for any improved action involving terms in which il) and ip are separated by more than a single link, there is no 
longer a local conserved baryon number current bilinear ji_i{x) such that ^ (in (a;) —3iJ.(x~ji)) — for non-zero lattice spacing. 
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by replacing d'^M/dfi^ with dM/dm = 1 and d^M/dm" = for n > 2. In the case of the reweighting for m, we 
compute the fermionic reweighting factor up to second order, and the chiral condensate up to first order, i.e. 

Indet Af (to) - Indet Af (toq) = trAf~^(TO - toq) - tr(A/^^Af"^)(TO - mof/2 + 0[{m ~ toq)^], (22) 

= {N^Nt)-^aNf[tTM-^ ~ tr{M-'^M~^){m ~ toq)] + 0[{m - toq)^]. (23) 

Hence, the error of the Polyakov loop susceptibility is 0[(to — toq)'^] and that of the chiral susceptibility 0[(to — toq)^]. 
Figures |l| and || show xl and Figs. || and || show for different to as functions of /3 for simulation masses toq = 0.1 
and 0.2. These figures show that the peak position moves to smaller /3 as to decreases as expected. Moreover we 
find that as to decreases the peak height becomes lower for xl and higher for Xvii/-- These behaviours are consistent 
since the Polyakov loop is an exact order parameter only in the limit to — > c», while the chiral condensate is an 
order parameter only for m — > 0. The phase transition is known to be a crossover for 2 flavor QCD with to > 0. 
We calculate the slope of the transition point (Mc/dm assuming that Pcifn) is defined by the peak position of the 
susceptibility whenever a clear peak is obtained^. Figures and |^ show /3c (w) for each toq. We fit the data by a 
power series expansion about toq, i.e., f3c{rn) — fSdnio) + J2n=i Cti(™ — wq)", with fit range |(to — too)/too| < 0.05 
or 0.1. The results are presented in Table ||. We find a linear fit to be adequate with the dependence on choice of 
A^fit less than 3%; the discrepancy from the choice of fit range is less than 10%. Both uncertainties lie well within 
the statistical error. We denote the fitted line for N^t = 1 and |(to — too)/too| < 0.1 by a dashed line. In Fig. ^ 
we compare the predicted variation of f3c{m) with previously existing data |14|. Filled symbols are the results of the 
current study. The short lines denote the upper and lower bounds on the slope (3'^. From this figure, we find that 
reweighting yields results which are quite consistent with those of direct simulation, and hence infer that reweighting 
the fermion action using the technique we have outlined works well. 

V. REWEIGHTING FOR CHEMICAL POTENTIAL 
A. Chemical potential dependence of the transition temperature 

Next we turn our attention to reweighting with respect to /i, with the Taylor expansion made about the simulation 
point /i = 0. First we calculate the derivatives of the transition point with respect to /i in the region of small /j, 
relevant to RHIC. In Figs. H, |l^, ^ and [ij, we plot xl and at to = 0.1 and 0.2 for various fi. As outlined in 
section |l[ we compute consistently up to 0{ii^) and expect the results to contain errors at 0{p^). Strictly speaking, 
the 0(/i^ term does not vanish for L since it is complex (see Sec. ||). However, we expect that xl and x^,/, yield the 
same /3c (see below) with error 0(/i*). The figures show that the position of the susceptibility peak moves lower as 
^ increases, which means that the critical temperature becomes lower as /i increases. As we obtained well-localised 
peaks for xl at to = 0.2 and x^bTb ^ — ^-^ ^'^^ 0.2, we use these peak positions to determine the transition point 
Pc as functions of /i^ in Figs. |l3|, and ^ Note that because the Polyakov loop is interpreted as an external quark 
current running in the positive time direction, positive and negative /i give different contributions to both L and xl, 
and we display both cases. Figs. |l^, |lj and |l5| also display the value oi ^ = O.lTc relevant for RHIC. The shift 
/3c(/i) — /3c(0) is found to be small at this point. 

Because the first derivative is expected to be zero as discussed above, we fit the (3c data by a straight line in /i^, fixing 
/3c at /i — 0, in ranges < 0.008(0.014) for m = 0.1(0.2) respectively, in which the phase problem is not serious (see 
Sec. |VC| below). Wc obtain d^/3c/d^^ — —1.20(44) and —1.02(56) at to = 0.1 and 0.2 from the chiral susceptibility 
and d^c/d/i^ — —1.01(23) at m = 0.2 from the Polyakov loop susceptibility. Dot-dashed lines in Figs. ^ |l^ and 
[l5| are the fitted lines. To investigate the fit range dependence and the fitting function dependence, we also tried the 



range /^^ < 0.005(0.01) for to = 0.1(0.2), and using a quadratic fit function. Table HI summarises the results. We 
may conclude that |d^/3c/d/i^| w 1.1 with 30 — 50% error, and any quark mass dependence of d^/3c/d/j,^ is not visible 
within the accuracy of our calculation. 

Of course, it is desirable to translate these observations into physical units. The second derivative of Tc can be 
estimated by 

d^Tc 1 d^/3c/^^d^\ (^^^ 



d^il N^Tc d^i^ / V da 



^Because the peak width of X-L is too wide for the smaller mass m — 0.1, we do not determine the pseudocritical point for L 
in this case. 
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where a is the lattice spacing. The beta-function may be obtained from the string tension data in Ref. [|^. We 
compute it by differentiating the interpolation function of the string tension with an ansatz : 

V^iP) = i?(/3)(l + C2a\p) + C4a*iP))/co, (25) 

where R{P) is the usual 2-loop scahng function, a = R{f3)/R{j3) and /3 = 3.70. Co,C2 and C4 are fit parameters 
with Co = 0.0570(35), C2 = 0.669(208) and C4 = -0.0822(1088) at m = 0.1. We get a-i(da/d/3) = -2.08(43) at 
(/3,to) = (3.65,0.1). We then find Tc{d^Tc/d^i\) w -0.14 at to = 0.1. We sketch the phase transition line with 50% 
error in Fig. |l^ assuming Tc — 170MeV. In the figure we also indicate the line /iq/T = 0.4, corresponding roughly to 
the range over which the fits to the leading order behaviour of T'c(/Lt) shown in Figs. |l^ ~|l^ are made. Of course, one 
has to expect that higher order terms in the expansion become relevant for ^/T — 0{1). To quantify this we will have 
to analyze higher order contributions in the expansion in the future. To indicate the present systematic uncertainty 
in the transition line for larger /i/T we show this region as a dotted line in Fig. We stress that the errors shown 
are statistical only and reflect the uncertainty of the coefficient of the Oip^) term in the expansion of Tc{fJ,). On the 
assumption that the transition line is parabolic all the way down to T = 0, then this curvature is too small to be 
consistent with the phenomenological expectation that at T = a transition from hadronic to quark matter occurs 
for /ic some 50 - 200 MeV greater than the onset of nuclear matter at ~ mAr/3 ~ 307MeV |lj]. This tendency 
is also supported by the result of Fodor and Katz [Q, and hints at contributions from higher-order derivatives, or 
even non-analytic behaviour, at larger values of /i. Despite the large errors we can see that our result gives us useful 
information about the phase diagram, at least for small /i, because the first derivative is zero. 

Another point worth noting is the screening effect of dynamical anti-quarks at ^ < 0. Negative chemical potential 
induces the dynamical generation of anti-quarks, which in contrast to quarks can completely screen an external color 
triplet current. Thus the free energy of a single quark is reduced, especially in the confinement phase, and the 
singularity at the phase transition point is weakened due to the reduction in the range of current-current interactions. 
This effect can be seen in Figs. ^, ^0|, |l^ and |l^, where we denote the Polyakov loop and its susceptibility at /i < 
by dot-dot-dash and dot-dash-dash lines. We see that L at /i < is larger than that at /i > 0, which means that the 
free energy at /i < is smaller. Moreover, as seen in Fig. |o[ the peak height of xl becomes smaller for /i < 0, while 
the position of the pseudocritical line in Fig. |l^ is almost the same between positive and negative /z. The screening 
effect only seems to make the phase transition singularity weaker without shifting the pseudocritical line. Because 
the only source of asymmetry between fi and — /i is due to the correlation between the imaginary parts of the fermion 
determinant and L, these imaginary contributions help to decrease the susceptibility at fi < 0. In this way, we can 
see that the explicit breaking of time reversal symmetry by exchange of fi with —fi helps to highlight an interesting 
feature of dynamical quarks in full QCD. 

Finally, if instead we were to impose an isovector chemical potential fij having opposite sign for u and d quarks 
1^,^, then the quark determinant would become real and positive, enabling simulations using standard Monte-Carlo 
methods |20|. This motivates a comparison between systems with the usual isoscalar chemical potential /i and the 
isovector chemical potential fij. In the framework of the Taylor expansion, terms even in fi are identical for both u 
and d quarks, but odd terms cancel for the case /x/ ^ 0, meaning that terms proportional to Oi,TZi should be set 



to zero in eqn. (14). We analyzed the transition point (3c{fJ-i) for m = 0.2; the results are shown in Fig. for Pc 
determined by xl and Fig. |2^ for that by x^^. The solid line shows (3c as a function of ^ij, the dashed line Pdp). 
The second derivative of (3c with respect to /x/ is found to be —0.96(19) for xl and —0.93(52) for Xtft^,- Dot-dashed 
lines in Figs p^and[2C| show the fits. Within errors there appears to be no significant difference between isovector and 
isoscalar chemical potentials for small /z. A similar analysis for at m = 0.1 is shown in Fig. here the second 
derivative of (3c is —0.71(16), which is smaller than the isoscalar case. However this result is also smaller than that 
obtained at m = 0.2, which is physically unacceptable since the second derivative should approach zero as m — s- 00. 
Hence the difference between fij and /i at to = 0.1 is most likely due to statistical error. 

The terms we have dropped are associated with fluctuations in the phase of det M , which are small in the region of 



small /i, as will be demonstrated in Sec. V C below. This is perhaps not unexpected on physical grounds - increasing 
/i/ is predicted to induce the onset of matter in the form of a pion condensate at a critical /i/o ~ mps/2 [ p^ , and 
indeed evidence for this scenario in the form of a negative curvature for mps{fJ-i) in the low T phase is reported in 
pof . However, even for to = 0.1 on this lattice this scale is roughly 0.92^/a ~ 390MeV [Q, which is a little larger 
than the isoscalar onset threshold Ho ^ WAr/3. The curvature with repect to /i/ should dominate as the chiral limit is 
approached and pion and nucleon mass scales become separate. If this turns out to be the case, then it is interesting 
to note that phase correlations between observable and measure actually decrease the physical effect of raising /x; this 
has also been observed in simulations of two-color QCD with a single flavor of staggered adjoint quark [^, in which 
including the sign of the fermion determinant has the effect of postponing the onset transition. 
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B. Quark number susceptibility and equation of state at ^ 7^ 



The energy density e and pressure p at the critical point are interesting quantities for heavy- ion colhsion experiments. 
In this section, we discuss the ^-dependence of the equation of state which describes them. If we employ the integral 
method based on the homogeneity of the system |^^, we obtain p = (T/V)hiZ; derivatives of p with respect to fi 
are then related to the quark number density nq (via a Maxwell relation) and the singlet quark number susceptibility 
Xs = driq/dfiq [|: 



1 d\nZ 
1 a^lnZ 



rp4 



Here Uq, xs and also the nonsingiet susceptibility xns are given in physical units by 



riq 


aNiN^ 


T3 


N! 


XS 


aNfNt 


y2 





I 1 dM\ 



w"" 



[aNifNt 



Xns 
T2 



aNfNt 

~Nr~ 



tr M 



tr 



0/1 o^i J 



(26) 
(27) 

(28) 

(29) 
(30) 



The quark number density is zero at /i = so once again the leading correction is 0{^^). The susceptibilities xs^i^j 
and xnsQ.^ are plotted in Figs. |2^ and ^ for m = 0.1 and 0.2. Because xs^^ = 0.0433(3) and 0.0306(2) for m = 0.1 
and 0.2 at (3c in Table |lT| (^/^V') , we obtain T^d^{p/T'^)/dfil = 0.693(5) (m = 0.1) and 0.490(4) (m = 0.2) at P^- The 
discrepancy of p/T^ at the interesting point for RHIC, fJ,q/Tc ^ 0.1, from its value at /i = is about 0.0035(0.0024) 
for TO = 0.1(0.2); since p/T^ « 0.27 at /3c for (TO,/i) = (0.1,0) ||l|l this is a 1% effect, and hence quite small. We 
can also obtain estimates of the quark number density via Uqa^ ~ //qoxsa^, with results Uq/T^ = 0.693(5)/iq/r and 
0.490(4)/Xq/T for m = 0.1 and 0.2 which assuming T ~ 170MeV translates into roughly 9% and 6% of nuclear matter 
density at the RHIC point. Clearly these values will need careful extrapolation to the chiral limit before a meaningful 
comparison with experiment can be made. 

Moreover the energy density e can also be estimated via the equation for the conformal anomaly: 



3p 



rp4 



1 d\nZ 
a 



1 



da 



dm 9 In Z 
da dm 



da dp 

Here we estimate e in the chiral limit, where adm/da can be neglected. We find 

e~3p 1 



rp4 



VT^ dl3 



dlnZ f I da 
adp 



with second derivative 



dHie^3p)/T^) 
d^il 



J_9xs 

't^ dp 



Ida 
adp 



(31) 



(32) 



(33) 



Because the quark mass dependence of the equation of state seems to be small in Ref. |23| , we estimate the derivative 
using the value of xs at m = 0.1 and 0.2. Using the formula, d{0)/dp = (0{~dS/dP)) - (0)(~dS/dP), we obtain 
d{xsO'^)/9P = 1.11(5) and 0.82(4) at Pc for to — 0.1 and 0.2. Then the second derivative of e — 3p is estimat ed to 
be T^d^ie 



Zp)|T^)/d^Jii = 8.5(1. 



Finally, we obtain T^d^{e/T^)/dni^ = 10.6(1. 



at TO = 0.1, where we use the same value of the beta- function as in Sec. V_A. 

The discrepancy of e/T^ at the RHIC point from ^ = is about 0.05. 
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Once again, because e/T" « 6 at /3c for (m,/i) = (0.1,0) ||, this is a 1% effect, suggesting that the /iq-dependence 
of the equation of state is small in the regime of relevance for RHIC. 

Next we discuss the relation between the equation of state and the phase transition line. It is of great interest to 
investigate whether the values of the pressure p(Tc(/iq), fiq) and energy density e{Tc{fJ.q), fj-q) along the transition line 
are constant or not. To this end, consider the line of constant pressure in the (T, fiq) plane, i.e.. 



dp 

dU4 



T 



d{p/T^ 
dT 



4p 



AT 



T 



d{p/T^) 



(34) 



together with a similar relation for Ae, and compare it with the phase transition line. The slope of the constant 
pressure line is then given by 



The derivative d{p/T^)/dT can be calculated by 



dr ^ d{p/T^) I (diprn , 4p 

d(A^?) ^{^,D / \ dT 



. d{p/T^) 
dT 



Iday^ d{p/T^) 



adf3 



d(3 



1 da 



I 1 dSa 



' a 9/3 j ^'\\ NfNt dp ] \ NlNt d/3 )J ' 



/ 1 dS, 



(35) 



(36) 



where means the expectation value evaluated at T = for normalization. Using the da ta o f Ref. [ p^ , 

p/T^ = 0.27(5), d(p/T^)/dl3 = 4.5(9) at T^ for m ^ 0.1, together with the beta-function in Sec. we obtain 



T{d{p/T'^)/dT)\T=T, = 2.2(6) for m = 0.1. Noting also that dip/T'^)/d{pl) = {l/2){d^{p/T'^)/dnl) = 0.347(3)/r2, 
we find that the slope of the constant pressure line emerging from the critical point on the T-axis is T(dT/d(/iq)) = 
—0.107(22). A similar argument using the data of [|4j gives the slope of the constant energy density hne T(dT/d(^q)) = 
-0.087(23). Because the slope of the transition line in terms of /i^ is Tc{dTc/d{fil)) = {l/2)Tc{(PTc/dnl) w -0.07(3), 
we deduce that the variations of p and e along the phase transition line are given by 



p(T,{^lq),^lq)-p{TM,0)=^4TH0) x 0.12(11) ; e{TM,fiq) - e{TM,0) ^ ^^lT^{0) x 1.0(2.2), 



(37) 



the dominant source of uncertainty in each case being the location of the phase transition line itself. Within our 
errors, therefore, both pressure and energy density appear constant along the phase transition line. 



C. The phase of the determinant at p / 



Finally we discuss the region of applicability of generic reweighting approaches. If the reweighting factor in eqn. (|i|) 
changes sign frequently due to the complex phase of the quark determinant, then both numerator and denominator 
of (|l|) become vanishingly small in the thermodynamic limit, typically behaving ~ e"^""" with the lattice size A'sitc = 
N^Nf. This makes control of statistical errors in the calculation of the expectation value very difficult. Of course, 
arg(det M) starts at zero at /i = but grows as /x increases. It is important to establish at which value of /i the sign 
problem becomes severe. 

As discussed in section ||, the phase can be expressed using the odd terms of the Taylor expansion of Indet M. If 
we write det A/ = I det A/|e'^, then 



9 = aNf Im 



SlndetM fi^d^lndetM 
^ d^ 3! 



(38) 



For small /i, the first term Q;iVtImtr[Af ^^(9Af3/i)]/i is dominant. Now, because {NlNt)~^ir[M^^{dM / dfj)] is 
the quark number density, its expectation value must be real and in fact vanishes at /i = 0. Although the 
average of the phase is zero, its fluctuations remain important. We investigated the standard deviation of 



{N'l Nt)~^ln\iv[M'^ {dM / d^j)] and present the results in Table |\]. We find values of about 2.2 x 10"^ at /3c(m = 0.1) 
and 1.6 x 10"'^ at /3c(™ = 0.2). The standard deviation of the leading term of ( ^ ) therefore has a magnitude of about 
18/i for m — 0.1 and 13/i for m = 0.2 in the vicinity of the transition. Consequently the phase problem appears from 
ji ^ 0.09(0.12), i.e., ^q/Tc ~ 0.4(0.5) for m — 0.1(0.2), since the phase problem arises if the phase fluctuation becomes 
of 0(1). We notice that the value of /i for which the phase fluctuations become significant decreases as either m or 
(} decreases. Roughly speaking, the numerator and denominator of (|^) decrease in proportion to the average of the 
phase factor (Re(e**)). We show this factor for various /3 and m in Fig. p3, where it is clear that the average becomes 
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small around the values of qu oted abo ve. T he phase fluctuations at the RHIC point /iq = O.lTc, however, are small 
enough for the analysis of sees. V A and VB to be applicable. 

We should also note that the fluctuation of the phase depends on the lattice size iVgitc, and on the number of the noise 
vectors N^- From general arguments, the phase of the reweighting factor is expected to decrease as (e**) oc e~^="°, 
implying that the applicable region of reweighting becomes narrower as the lattice size grows. By contrast, the value 
of Im tr\AI~^{dM/dfi)] calculated on each conflguration also contains an error due to the finite number of noise 



vectors (see eqn.( 
as seen in Table 



of the Appendix); for Nn = 10 this error is not small compared to the standard deviation, 
The phase fiuctuation discussed above includes this error due to finite N,^, and we suspect 
that the true fiuctuation becomes smaller as TVn increases. To confirm this, we reanalyze the standard deviation 

^({Imtr[M-i(aM/5Ai)]}2) - {lmtT[M-^^M/^^l)]f by treating the calculation of {{lintT[M-\dM/dn)]}^) more 
carefully. Since the noise sets must be independent, we subtract the contributions from using the same noise vector 
for each factor. Details are given in the Appendix. The results are quoted in the STD(Imp.) column of Table IV 
and are found to be significantly smaller. Because they might be closer to the = oo limit, they suggest that the 
standard deviation for larger A^n is much smaller, which means that the region of applicability becomes wider as A'n 
increases. 



VI. CONCLUSIONS 



In this paper we have proposed a new method based on a Taylor expansion in chemical potential fi to investigate 
the thermodynamic properties of QCD with /i ^ 0. By computing the chiral susceptibility and the Polyakov loop 
susceptibility for 2 flavors of p-4 improved staggered fermions, we have been able to estimate the dependence of /3c, 
and hence the critical temperature Tc, on fi on moderately large volumes, thus reinforcing the recent advance of 
lattice QCD into the interior of the (/Ltq,T) plane [Q. We have also been able to quantify the effect of a non-zero 
chemical potential on the equation of state. Although we have focussed on critical observables in order to fix physical 
scales, the method can be applied in a small range of /i at arbitrary (3, although the radius of convergence is expected 
to decrease as T — s- since in this limit all /^-dependence should vanish for /Xq < fig, making the behaviour about 
the origin non-analytic. The method is also applicable to a range of physical observables We find that Tc 

decreases as /i increases, but this appears to depend only weakly on quark mass, an effect also observed in studies of 
the equation of state p{T) |^^. Our results are in broad agreement with estimates based on exact reweighting 0| and 
suggest that the discrepancy of Pc from its value ai fj, — is small in the interesting region for heavy-ion collisions. 
Moreover we have observed evidence that when a negative chemical potential is imposed, the generation of dynamical 
anti-quarks and the consequent screening of an external color triplet current is enhanced. 

An unresolved issue is the method's limitations. We have been able to estimate the complex phase of the fermion 
determinant for a 16^ x 4 lattice and found that the sign problem is not serious in the range iiq/T^ < 0.4-0.5 for 
m = 0.1-0.2, covered by this study. It is not yet clear to us to what extent the radius of convergence of the Taylor 
expansion is linked to the fluctuations of arg(detM). An optimist might hope that the method can yield accurate 
thermodynamic information all the way out to the critical endpoint where the quark/hadron phase transition changes 
from second to first order; moreover, since individual terms in the expansion are expectation values of local operators, 
the method should be applicable on arbitrarily large volumes, particularly if larger numbers of stochastic noise 
vectors than we have used here are employed. A pessimist might worry that phase fluctuations should make calculation 
of higher order terms impracticable long before the radius of convergence is reached, particularly as the chiral limit 
is approached since in this case the correlations between arg(det M) and Im(C') should discriminate between the 
different physics associated with isoscalar and isovector chemical potentials. More work is needed before we can say 
which is more realistic. 

After this work was submitted we learned of a paper which studies the phase transition line by analytical continuation 
of results obtained by simulation with imaginary fi p5| . The results are in reasonable agreement with ours. 
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APPENDIX: REMARK ON THE NOISE METHOD 



The calculation of an operator such as (trv4)^, where A is a matrix, using the noise method has to be treated 
carefully. Because the random noise vectors should be independent for each calculation of tr^, 



(trA)2 



1 ^" 1 ^" 1 

]^ E E - ^i™^ jV„(iV„_l) E 

" a=l " 6=1 " ' a^b 



This equation can rewritten as 



(trAY = lim 



where e{A) is the error due to finite A^n: 



(39) 



(40) 



(41) 



The error decreases as {Nn — 1)~^ as A^n mcreases, but can be significant for small N^- Moreover, e'^{A) is negligible 
for an operator which always has the same sign such as trM~^; in this case its contribution is about 0.001% for 
((trM^^)^) with Nn = 10. However, for an operator which changes sign frequently such as tr[Af~^(9Af/9/i)], the 
effect of the additional term is important. We calculate the quark number susceptibility and the value of 'STD(Imp.)' 
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in Table IV taking this additional term into account. The difference between 'STD' and 'STD(Imp.)' in Table IV is 
the contribution from the additional term. 

Next, we construct the reweighting method based on Taylor expansion, eqn.(^, explicitly up to second order using 
the noise method. Assuming O is a bosonic operator, we can rewrite the numerator of eqn.(y): 



J 



Otr M 



1 dM 1 dM\ 
OtriM'^—M-^—) 



(42) 




0( r/tM-i^r;^ 



- ( O UtM 



dM 
djjL 



■M- 



9A/ 
9/i 



(43) 



= ( Oexp <( /iaiVt^77tM-i^77^ 



(44) 



where (• • •) denotes the average over the noise vectors. The denominator of eqn.(|2|) is given by the same expression 
with O — 1. In each case a term proportional to appears. In Fig. ^ we estimate the effect of this term by 
subtracting this term from original one. The difference for xl by the subtraction is found to be quite small, e.g., that 
is less than 1% at m = 0.2 and /i < 0.1. The result suggests the contribution from the term of is small for xl 
although the value of e[M~^{dM / dy)Y itself is_not small. 

For the case of a fermionic operator such as many such additional terms appear in the reweighting formula. In 
this study, we neglect the effect from further additional terms, since Fig. |5] suggests that the effect is small for the 
determination of /3c. 



TABLE I. Simulation point (m, 13) and number of configurations A^conf for mass-reweighting and /i-reweighting. 

P A''conf(maSs) iVconf(/i) 



0.1 


3.640 


38000 


20000 




3.645 


15000 






3.650 


58000 


38000 




3.655 


16800 






3.660 


55000 


40000 




3.665 


7800 






3.670 


30000 


30000 


0.2 


3.740 


5000 






3.750 


30000 


20000 




3.755 


15000 






3.760 


52000 


34000 




3.770 


48000 


32000 




3.780 


5000 
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TABLE II. Quark mass 


dependence of 


transition point determined by 


L and {ipip). The fitting 


function is 


/3c = /3c(mo) + Er=iC.(m- 


mo)" The truncation error is contained in C2 from iptp. 




rno 


/3c (mo) 


Cl 


C2 


fit range 


AT 

JVflt 


^pij) 0.1 


3.6492(22) 


1.05(14) 


- 


-0.01 < m - mo < 0.01 


1 




3.6492(22) 


1.03(13) 


[-9.(14)] 


-0.01 < m - mo < 0.01 


2 




3.6492(22) 


1.07(19) 




-0.005 < m - mo < 0.005 


1 




3.6492(22) 


1.07(19) 


[-17.(26)] 


-0.005 < m - mo < 0.005 


2 


0.2 


3.7617(36) 


0.896(90) 




-0.02 < m - mo < 0.02 


1 




3.7617(36) 


0.894(89) 


[5.(13)] 


-0.02 < m - mo < 0.02 


2 




3.7617(36) 


0.970(168) 




-0.01 < m - mo < 0.01 


1 




3.7617(36) 


0.999(180) 


[18.(39)] 


-0.01 <m-mo <0.01 


2 


Polyakov 0.2 


3.7639(19) 


0.838(64) 




-0.02 < m - mo < 0.02 


1 




3.7639(19) 


0.835(63) 


-2.7(4.5) 


-0.02 < m - mo < 0.02 


2 




3.7639(19) 


0.883(106) 




-0.01 < m - mo < 0.01 


1 




3.7639(19) 


0.885(106) 


-4.7(10.0) 


-0.01 < m - mo < 0.01 


2 


TABLE III. /3c and its second derivative with respect to /i. We fit the data with the function /3c(m) = /3c(0) + y^^Ii, c„ij,^", 


where d^/3c/d/u^ = 2ci. 












m 


/3c 






fit range iVfit 




0.1 


3.6497(16) 




-1.20(44) 


< < 0.008 1 






3.6497(16) 




-1.19(54) 


< /i^ < 0.005 1 






3.6497(16) 




-1.21(79) 


< At^ < 0.008 2 




0.2 


3.7641(37) 




-1.02(56) 


< < 0.014 1 






3.7641(37) 




-1.10(68) 


< //^ < 0.010 1 






:!.7()11(:!7) 




-i.:il(i():!) 


< //- < O.Dii 2 




Polyakov 0.2 


3.7651(16) 




-1.01(23) 


< //"^ < 0.014 1 






3.7651(16) 




-1.07(24) 


< At^ < 0.010 1 






3.7651(16) 




-1.21(31) 


0<n^< 0.014 2 




TABLE IV. Average of (Ii 


m tI[{^M/^^J,)M~'^]), average < 


of its error for each configuration ((e)), standard deviation (STD) 


and improved standard deviation (STD(Imp.)). 










in J 


(Im tr[((/i\7/ 


OlOM ']) 


W) 


STD 


STl)(lmp.) 


0.1 3.64 


-1.15 X 


10-* 


0.00199 


0.00233 


0.00110 


3.65 


1.02 X 10"® 


0.00194 


0.00223 


0.00099 


3.66 


-3.06 X 10"® 


0.00189 


0.00212 


0.00085 


3.67 


-1.40 X 10"® 


0.00185 


0.00206 


0.00077 


0.2 3.75 


1.03 X lO"'' 


0.00141 


0.00168 


0.00085 


3.76 


0.93 X 10"® 


0.00140 


0.00161 


0.00072 


3.77 


-4.17 X 10"® 


0.00138 


0.00155 


0.00061 
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FIG. 1. Quark mass dependence of xl as a function of B at. mo = 0.1. 




FIG. 2. Quark mass dependence of xl as a function of /3 at mo = 0.2. 
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FIG. 5. I3c{m) determined by xl around mo = 0.2. 
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FIG. 6. I3c{m) determined by x^iV" around mo =0.1. 
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FIG. 7. I3c{m) determined by Xv;v> around mo = 0.2. 
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FIG. 8. (3c(rn) determined by Xi/.^ in comparison with previous results. 
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FIG. 9. xl{I3) at m = 0.1 for various /u. 
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FIG. 10. Xl{I3) at m = 0.2 for various jj,. 
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FIG. 13. Phase transition point /3c (m) determined by at m = 0.2. 
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FIG. 14. Phase transition point /3c(/i) determined by Xi,^ji at m = 0.1. 
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FIG. 15. Phase transition point Pc{^) determined by Xii^p a.t m — 0.2. 
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FIG. 16. Sketch of the phase diagram, as estimated using our value of the curvature of /3c (m = 0). The errors shown 
are statistical only and reflect the uncertainty of the coefflcient of the 0(/^^) term in the expansion of Tc(p). Dotted line is 
^/T — 0.4. The diamond symbol is the end point of the first order phase transition obtained by Fodor and Katz [M]. 
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FIG. 19. Difference between fj, and /// for /3c determined by Xi at m = 0.2. 
3.768 I > 



3.766 



3.764 



o 
CO. 



3.762 



3.760 



3.758 



3.756 



T 



T 




o.ood' 

RHIC 



Iso-scalar 

Iso-vector 

fitted line 



0.005 



0.010 



FIG. 20. Difference between /j, and /uj for /3c determined by XiAV at m = 0.2. 
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FIG. 21. Difference between /j, and fii for f3c determined by Xi/,v at m = 0.1. 
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FIG. 22. Quaxk number susceptibilities xs and xns at m = 0.1. Ax = Xs — Xns- 
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FIG. 25. Effect from the term of 0{£^) on xl at m = 0.2. Solid lines are the same as Fig. |l^ obtained including the O(e^) 
term, and dashed lines are calculated without it. 
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